/******************************************************************************
* FILE: mm.c
* DESCRIPTION:  
*   In this template code, the master task distributes a matrix multiply
*   operation to numtasks-1 worker tasks.
*   NOTE1:  C and Fortran versions of this code differ because of the way
*     arrays are stored/passed.  C arrays are row-major order but Fortran
*     arrays are column-major order.
* AUTHOR for MPL version: Ros Leibensperger / Blaise Barney
* LAST MPL REVISED: 09/14/93 for latest API changes.  Blaise Barney
* CONVERTED TO MPI: 11/12/94 by Xianneng Shen 
* MODIFIED: 4/17/09 by Joe Antoon
******************************************************************************/

#include <stdio.h>
#include "mpi.h"

#define NRA 500
#define NCA 500
#define NCB 500

int main( int argc, char *argv[] )
{
    MPI_Status status;
    int rank, size, workers;
    int src, dst, rows;
    int i, j, k;
    double start, sending, end;

    double a[NRA][NCA];
    double b[NCA][NCB];
    double c[NRA][NCB];

    int printResult = 0;

    MPI_Init( &argc, &argv );
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);
    workers = size-1;

    start = MPI_Wtime();

    if( rank == 0 )
    {
        // Allocate array
        for( i=0; i<NRA; i++ )
            for( j=0; j<NCA; j++ )
                a[i][j]= i+j;
        for( i=0; i<NCA; i++ )
            for( j=0; j<NCB; j++ )
                b[i][j]= i*j;
    }

    sending = MPI_Wtime();

    MPI_Bcast( b, NCA*NCB, MPI_DOUBLE, 0, MPI_COMM_WORLD );

    if ( rank == 0 )
    {
        // Send matrix
        int offset = 0;
        for( dst=1; dst<size; dst++ )
        {
            rows = NRA / workers;
            if( dst <= NRA % workers )
                rows++;

            MPI_Send( &rows, 1, MPI_INT, dst, 0, MPI_COMM_WORLD );
            MPI_Send( &a[offset][0], rows*NCA, MPI_DOUBLE, dst, 0, MPI_COMM_WORLD);
            offset += rows;
        }

        // Receive results
        offset = 0;
        for( src=1; src<size; src++ )
        {
            rows = NRA / workers;
            if ( src <= NRA % workers )
                rows++;

            MPI_Recv( &c[offset][0], rows*NCB, MPI_DOUBLE, src, 0, MPI_COMM_WORLD, &status );
            offset += rows;
        }

        end = MPI_Wtime();

        // Print results
        if ( printResult )
        {
            int i,j;
            for ( i=0; i<NRA; i++ )
            {
                for ( j=0; j<NCB; j++ )
                    printf("%.0f,", c[i][j]);
                printf("\n");
            }
        }

        printf( "\nFILLING: %fs\n", sending-start );
        printf( "CALCULATING: %fs\n\n", end-sending );
    }
    else
    {
        // Receive
        MPI_Recv( &rows, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &status );
        MPI_Recv( a, rows*NCA, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status );

        // Calculate
        for( k=0; k<NCB; k++ )
        {
            for( i=0; i<rows; i++ )
            {
                c[i][k] = 0.0;
                for( j=0; j<NCA; j++ )
                    c[i][k] += a[i][j] * b[j][k];
            }
        }

        // Send results
        MPI_Send( c, rows*NCB, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD );
    }

    MPI_Finalize();
}

